


library(mvtnorm)
library(arm)
library(BRugs)
library(R2OpenBUGS)
library(coda)
library(car)
library(foreign)
library(vioplot)


# Change this to your working directory.  Note this path makes Dropbox work at least on a windows machine
setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S6")


dataset<-read.dta(file.choose()) # read in a Stata file "ireds_dataset_wide_v12.dta"



# Data preprossing for the Congdon model 


# create the matrix of lab adjacencies

LabAdj=matrix(data=0, nrow=length(dataset$labid), ncol=length(dataset$labid))
for (i in 1:length(dataset$labid)) {
  for (j in 1:length(dataset$labid)) {
    if (i!=j) { 
      if (dataset$labid[i]==dataset$labid[j]) {
        LabAdj[i,j]<-1
      }}}}



#### Create Vector of labmate mappings to read into Congdon programs: 

LabMap<-list(map=NULL)
for (i in 1:length(dataset$labid)) {
  for (j in 1:length(dataset$labid)) {
    if (LabAdj[i,j]==1) LabMap<-list(map=c(LabMap$map, j))
  }}

LabC<-list(C=0)
nadj<-0
for (i in 1:length(dataset$labid)) {
  nadj<-sum(LabAdj[i,])+nadj
  LabC<-list(C=c(LabC$C, nadj))
}
map <- LabMap$map
C <- LabC$C

one<-rep(1, length(dataset$labid))
numLabNeigh<-sum(as.numeric(t(LabAdj%*%one)))
Cn <- numLabNeigh




## create data 


# set constants 
n_respondents <- length(dataset$labid)
n_responses <- 4

n_q6 <- 5

# create matrices of the pre and post outcomes.  You can include an arbitrary number of items/questions provided
# you have at least three to identify the latent variables.
S6.0 <- as.matrix(cbind(dataset$t0_labdyn_disagreevalid, dataset$t0_genrcr_discussdominate, dataset$t0_labdyn_discussrespect, dataset$t0_labdyn_discussunderstand, dataset$t0_genrcr_coworkerspractice))
S6.1 <- as.matrix(cbind(dataset$t1_labdyn_disagreevalid, dataset$t1_genrcr_discussdominate, dataset$t1_labdyn_discussrespect, dataset$t1_labdyn_discussunderstand, dataset$t1_genrcr_coworkerspractice))

# choose a covariate matrix:
X <- as.matrix(cbind(dataset$urm, dataset$male, dataset$pi_role, dataset$urmxcond, dataset$malexcond, dataset$pi_rolexcond))
#X <- as.matrix(cbind(dataset$lifesci, dataset$engineering, dataset$lifescixcond, dataset$engxcond))

# create variables that track (complete case) pretreatment means, which we use to improve convergence
s6x <- c(2.96, 3.28, 3.44, 3.17, 3.09)

condition <- as.numeric(dataset$condition)

cuts<-c(-1,0,1)




data<-list("S6.0", "S6.1", "n_q6", "s6x", "condition", "n_respondents", "n_responses", "C", "map", "Cn", "X") 
inits <- function() {
  list(S6.theta0=rnorm(n_respondents), S6.delta.theta=rnorm(n_respondents), 
       S6.alpha1=rnorm(1), S6.lambda.1=runif(n_q6),  
       S6.beta.1=rlnorm(n_q6), S6.beta.2=rlnorm(n_q6), 
       S6.lambda.0=as.matrix(rbind(cuts, cuts, cuts, cuts, cuts)),
       S6.beta.0=as.matrix(rbind(cuts, cuts, cuts, cuts, cuts)), S6.rho=0.1, 
       e6=rnorm(n_respondents), g=rnorm(length(X[1,]))
  )} 


# Output the data and inits files to the working directory for the analysis...
setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S6/wcovariates")
bugs.data(data, dir=getwd(), digits=5, data.file="data.txt")
bugs.data(inits(), dir=getwd(), digits=3, data.file="inits1.txt")
bugs.data(inits(), dir=getwd(), digits=3, data.file="inits2.txt")
bugs.data(inits(), dir=getwd(), digits=3, data.file="inits3.txt")

